Low-density neutron matter 
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The properties of low-density neutron matter are important for the understanding of neutron star 
crusts and the exterior of large neutron-rich nuclei. We examine various properties of dilute neutron 
matter using quantum Monte Carlo methods, with s- and p-wave terms in the interaction. Our 
results provide a smooth evolution of the equation of state and pairing gap from extremely small 
densities, where analytic expressions are available, up to the strongly interacting regime probed 
experimentally and described theoretically in cold atomic systems, where fcf « 0.5 fm" 1 and the 
pairing gap becomes of the order of magnitude of the Fermi energy. We also present results for 
the momentum distribution and pair distributions, displaying the same evolution from weak to 
strong coupling. Combined with previous quantum Monte Carlo and other calculations at moderate 
densities, these results provide strong constraints on the neutron matter equation of state up to 
saturation densities. 

PACS numbers: 21.65.-f, 03.75.Ss, 05.30.Fk, 26.60.-c 



I. INTRODUCTION 



The equation of state and the pairing gap of neutron 
matter at low densities are important for describing the 
properties of the inner crusts of neutron stars and pro- 
vide significant constraints for density-functional theories 
of large neutron-rich nuclei. Equation of state results 
at larger densities (p > 0.04 fm -3 ) have been used for 
some time to constrain Skyrme and other density func- 
tional approaches to large nuclei [HQ . More recently, the 
density-dependence of the 1 So gap in low-density neutron 
matter has been used to constrain Skyrme-Hartree-Fock- 
Bogoliubov treatments and especially their description 
of neutron-rich nuclei Q. At extremely low densities 
the equation of state and pairing gap can be expressed 
as analytically known functions of (fcpci), the product of 
the Fermi momentum and the neutron- neutron scattering 
length. Our results provide a smooth connection between 
these analytic results and previous calculations of neu- 
tron matter at larger densities, where the gaps become 
smaller and the superfiuid properties are less relevant to 
the equation of state. 

The properties of low-density neutron matter are par- 
ticularly important for describing the inner crust of a 
neutron star, which is composed of a lattice of neutron- 
rich nuclei along with a gas of neutrons and electrons. 
The neutron gas at low densities is expected to be su- 
perfiuid; the evolution of the equation of state and the 
pairing gap will impact the static and dynamic proper- 
ties of the inner crust of the neutron star. A cold neutron 
star will have a temperature from 10 6 K to 10 9 K (~ 0.1 
keV to 0.1 MeV); hence the low-density neutron gas is 
superfiuid because the critical temperature is expected to 
be larger, approximately 10 10 K (~ 1 MeV). The most 
basic aspects of a superfiuid gas embedded in a lattice 
of nuclei are described by zero-temperature infinite pure 
neutron matter. Corrections to this picture include gra- 



dient terms in the density induced by the ionic lattice. 
While these corrections are also important for density- 
functional theories of nuclei, we leave their study to fu- 
ture investigations. 

Superfluidity in neutron matter is often connected to 
cooling observations of neutron stars: the specific heat in 
a superfiuid is exponentially suppressed, a fact which is 
consistent with observations of cooling quiescent neutron 
stars [||. Furthermore, in the presence of a neutron 1 So 
gap, the neutron-neutron brcmsstrahlung reaction rate 
is also suppressed. Cooper-pair breaking/formation neu- 
trino emission processes that occur near the transition 
temperature are also relevant to the cooling of neutron 
stars during the crust's thermal relaxation |5|, @j. While 
many of these phenomena are not critically dependent 
upon the magnitude of the gap, recent heat-conduction 
mechanisms in magnetars require superfiuid phonons and 
their interaction with the lattice ions Q- These may be 
more sensitive to the magnitude of the gap. 

The neutron matter equation of state [8l4l 7| and pair- 
ing gap [l8l - l30j have been the subject of many studies 
over the years, with quite different results, particularly 
for the pairing gap. Now, however, cold-atom experi- 
ments can mimic the properties of dilute neutron mat- 
ter, giving nearly direct constraints on its properties. In 
cold atoms the interaction can be tuned through Fesh- 
bach resonances to produce a specific scattering length a, 
while the effective range r e between the atoms is nearly 
zero. In low-density neutron matter, on the other hand, 
the particle-particle interaction has a scattering length 
which is very large, ~ —18.5 fm, much larger than the 
typical separation between neutrons. The effective range 
is much smaller than the scattering length, r e ks 2.7 fm, 
so \r e /a\ w 0.15, but only at very low densities is the 
effective range much smaller than the interparticle spac- 
ing. We directly compare results in neutron matter and 
cold atoms to try to understand the impact of the effec- 
tive range theoretically; it may also be possible to use 
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narrow and wide resonances in cold atoms to study this 
experimentally. 31 

In a recent article [32| . we examined the similarities 
of cold atoms with neutron matter by calculating the 
T = equations of state and pairing gaps. The interac- 
tion used for the cold atoms had an infinitesimal range 
and the scattering length was varied to obtain results 
from fcpa = —1 to —10. For the case of neutron mat- 
ter, we took the s-wave part of the AV18 [33| interaction 
that fits s-wave nucleon-nucleon scattering very well at 
both low- and high-energies. In this work, we extend our 
approach to include p-wave interactions, examining their 
effects on the equation of state and supcrfluid pairing 
gap. Additional corrections due to higher partial waves 
and three-nucleon interactions are expected to be quite 
small in this density regime. We also calculate additional 
properties of neutron matter, including the quasiparti- 
cle excitation spectrum, the momentum distribution, and 
pair-distribution functions. 

All calculations are performed using quantum Monte 
Carlo (QMC) techniques (section lllip . including Varia- 
tional Monte Carlo (VMC) and Green's Function Monte 
Carlo (GFMC) methods. We compare our results to ana- 
lytic calculations at very low densities (section [II A\ , and 
to BCS calculations over the range of densities we con- 
sider (section III B[) . The BCS calculations are also used 
to try to understand and estimate the finite-size effects in 
the QMC simulations. Although the BCS results are not 
expected to be quantitatively accurate, they do provide 
a useful benchmark for comparisons and for understand- 
ing physical effects beyond the mean-field treatment of 
pairing. 



However, as was shown in 1961 by Gorkov and Melik- 
Barkhudarov,[35| the BCS result acquires a finite polar- 
ization correction even at weak coupling, yielding a re- 
duced pairing gap: 
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The polarization corrections reduce the mean-field BCS 
result by a factor of 1/ (Ae) 1 ^ 3 re 0.45. Interestingly, if one 
treats the polarization effects at the level of sophistica- 
tion used in the work of Gorkov and Melik-Barkhudarov, 
this factor changes with kpa (36[, though there is no a 
priori reason to expect such an approach to be valid at 
stronger coupling (kpa of order 1 or more). Calculating 
the pairing gap in this region has been an onerous task, 
as can be seen from the multitude of publications devoted 
to this subject in the past few decades. [l8l - l30j 



B. BCS in the continuum and in a box 

As the coupling strength increases, we expect the BCS 
mean-field theory to become more accurate. In the BCS- 
BEC transition studied in cold atoms, the BCS theory 
goes correctly to the two-body bound state equation in 
the deep BEC regime. Though we do not expect BCS re- 
sults to be precise, BCS theory provides a standard basis 
of comparison for our ab initio results and also allows 
us to analyze finite-size effects in the QMC simulations 
in a simple way. Within the BCS formalism the wave 
function is: 
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II. ANALYTIC RESULTS 

A. Weak coupling 

At extremely low densities (l^f a| << 1) the effective 
coupling between neutrons is weak and neutron matter 
properties can be calculated analytically. The ground- 
state energy of normal (i.e. non-supcrfluid) matter in 
this regime was calculated by Lee and Yang in 1957: (3~H 



where ut 



= 1. A variational minimization of the ex- 



pectation value of the Hamiltonian for an average particc 
number (or density) leads to the gap equation: 
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where the elementary quasi-particlc excitations of the 
system have energy: 
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where Efg is the energy of a free Fermi gas at the same 
density as the interacting gas. While this expression ig- 
nores the contributions of superfluidity, these are expo- 
nentially small in (1/kpa). In the next section we com- 
pare these results to QMC calculations for \kpa\ > 1. 

The pairing gap at weak coupling is also known ana- 
lytically. The mean-field BCS approach described below 
[Eq. [5] reduces in this limit to: 
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E(k) = v/£(k)2 + A (ky 
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and £(k) = e(k) — /i, where the chemical potential is \x 
and e(k) = is t ne single-particle energy of a particle 
with momentum k. The chemical potential is found by 
solving the gap equation together with the equation that 
provides the average particle number: 



i 



E(k) 



(7) 



When interested in the gap for neutron matter, it is 
customary to perform partial- wave expansions of the po- 
tential and the gap functions, as well as an angle-average 
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FIG. 1: (color online) BCS neutron-matter pairing gap A di- 
vided with the Fermi energy Ef, versus the Fermi momentum 
kF for AV18 (solid line) and a modified Poschl-Teller poten- 
tial (dashed line) tuned to have the same scattering length 
and effective range as AV18. At low density the two curves 
are identical for practical purposes. Also shown is the solu- 
tion of the BCS problem in a periodic box using the modified 
Poschl-Teller potential for 66 particles (dotted line). 



where vq and v are parameters which we tuned so that 
this potential reproduces the neutron-neutron scattering 
length a « —18.5 fm and effective range r e w 2.7 fm. 
The potential in Eq. (JTTJ clearly has no repulsive core, 
making it more amenable to a straightforward iterative 
solution. In Fig. Q] we show the results for these two 
potentials. For all the densities considered in this work, 
the results of solving Eqs. © and (|10p with these two 
potentials are virtually indistinguishable. For treatments 
beyond the mean field, though, more care must be taken. 
A simple purely attractive interaction with finite positive 
effective range will produce a collapse to a system size 
of the range of the potential. The repulsive core avoids 
this collapse in QMC calculations, but the details of the 
core interaction are not important at the low densities 
considered here. 

The modified Poschl-Teller potential can also be used 
in a calculation for finite average particle number. We 
have solved Eqs. (J5J and for (AT) from 20 to 200, in 
periodic boundary conditions in a cubic box of volume 
L 3 : 



k n = -j-{n x ,n v ,n z ) . (12) 



approximation. Thus, Eq. ([5]) takes the form: 

oo 

o 

where the potential matrix element is: 

oo 

v(k,k') = J dr r 2 jo(k'r)V(r)jo(kr) . 
o 
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Similarly, Eq. ([7J becomes: 
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These equations are one dimensional, and thus simple to 
treat numerically. The density equation can be decoupled 
from the gap equation only when A/ fx << 1; this is not 
the case for the density regime we are considering. 

We have solved Eq. © in tandem with Eq. (|10l) for 
the 1 Sq channel of the Argonne v!8 33j potential that 



contains a strong short-range repulsion. This calculation 
is greatly simplified if one uses the method described in 
Ref. [37j , thereby transforming the problem into a quasi- 
linear one. We have also solved Eq. ([5]) together with 
Eq. (fT0|) for a modified Poschl-Teller potential: 
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The solution to this problem for many values of (N) at 
a fixed density of kpa = —10 was given in Ref. [32j. 
There it was found that (N) = 66 is very close to the 
thermodynamic limit. We have performed similar calcu- 
lations for other densities, finding that they all exhibit 
the same trend. We have also performed similar compu- 
tations with generalized boundary conditions including 
separate phase shifts for spin up and down particles at 
the box boundary. [38| In Fig. Q] we show the results of 
solving the BCS gap equation Eq. (|5|) in a periodic box 
along with the particle-number conserving Eq. ([7J for 
(N) = 66. We do not expect this procedure to be suf- 
ficient for very weak coupling, \kpa\ < 1, as the pair 
size becomes larger than the simulation volume. A more 
detailed study is warranted to attempt to extract pair- 
ing gaps in this regime. In the rest of this work, when 
we mention BCS results in a box we will refer to the 
(N) = 66 case. 

We have also calculated, both in the continuum and 
in a periodic box, the momentum distribution, which in 
BCS is given by the following expression 
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as well as the energy of the quasi-particle excitations, 
which follows from Eq. (JS]). These are given in sections 
IllTEl and lllTFl 
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III. QUANTUM MONTE CARLO 
A. Hamiltonian 

The Hamiltonian for neutron matter at low densities 



is: 



N 



k=l 



(14) 



where N is the total number of particles. The neutron- 
neutron interaction is generally quite complicated, hav- 
ing onc-pion exchange at large distances, an intermediate 
range spin-dependent attraction by two-pion exchange, 
and a short-range repulsion. In the regime of interest, 
though, the scattering length and the effective range are 
most crucial to the physical properties of the system. 
For the purposes of our simulation, a short-range repul- 
sive core is also important so as to avoid a collapse to a 
higher-density state. 

In Ref. [23 we used the 1 So potential as the inter- 
action between all opposite-spin pairs. A perturbativc 
correction was added to correct for the fact that the 
S = I, Ms = pairs must be in a relative p-state (or 
higher) due to anti-symmetry. The p-wave interaction 
was neglected in those calculations. Here we improve this 
treatment by explicitly including p-wave interactions in 
the same-spin pairs, and perturbatively correcting the 
S = I, Ms = pairs to the correct p-wave interac- 
tion. We use the AV4 potential to determine the p-wave 
interactions. [39| . 

The AV4 interaction for neutrons can be simplified to 



Vi(r) = v c (r) + u ff (r)<ri • <x 2 , 



(15) 



which in the case of the S=0 singlet pairs gives 

v s (r) = v c (r) - 3v a (r) . (16) 

In this paper we add the contribution from spin 1 (triplet) 
pairs: 



v P (r) = v c (r) + v a (r) 



(17) 



The same-spin potential contribution is small even at the 
highest density considered. While still keeping the poten- 
tial of Eq. (jTB"]) in the propagator of our QMC method for 
the opposite-spin pairs, we have corrected perturbatively 
using the general case described in Eq. (|15[) . which can 
be written as (see, e.g., Ref. [40|): 

Vi(r) = v c (r) + v a (r)(-2P M - 1) (18) 

in terms of the Majorana exchange operator. 

B. Variational and Green's Function Monte Carlo 

We employ standard Variational and Green's function 
Monte Carlo methods to calculate the properties of dilute 



neutron matter. VMC calculations use Monte Carlo inte- 
gration to minimize the expectation value of the Hamil- 
tonian: 



{H)vMC = /dR|vMR)| 2 " 



(19) 



thereby optimizing the variational wave function \EV- 

Fixed-node GFMC simulations project out the lowest- 
energy eigenstate "to from a trial (variational) wave func- 
tion This they do by treating the Schrodinger equa- 
tion as a diffusion equation in imaginary time r and 
evolving the variational wave function up to large r. 
The ground state is evaluated from: 



*o = exp[—(H - E t )t]^ v 

= JJexp[-(F - E t )At]^ v , 



(20) 



evaluated as a branching random walk. The short-time 
propagator is usually taken as 

exp[-i? At] = cxp[-VAr/2] cxp[-TAr] exp[-VSr/2pl) 

which is accurate to order (At) 2 . For the lowest densities 
considered, we include the two-body propagator exactly: 



exp[-#Ar] = cxp[-TAr] 



cxp[— H 2 At] 
exp[— H At] ' 



(22) 



where the two-body Hamiltonian H2 includes the pair 
relative kinetic energy and the pair potential and Hq is 
the pair kinetic term only. At lowest order in (At) this is 
equivalent to Eq. 1211 However it includes multiple scat- 
tering for a pair and allows accurate calculations with 
larger time steps At. This is particularly important for 
very dilute systems where these multiple-scattering con- 
tributions of individual pairs dominate. 

The fixed node calculation gives a wave function ^q 
that is the lowest-energy state with the sames nodes (sur- 
face where ^ = 0) as the trial state "Jy. The resulting 
energy Eq is an upper bound to the true ground-state 
energy. The variational wave function "Jy has a Jastrow- 
BCS form (see next sub-section), and contains a variety 
of parameters, many of which affect the nodal surfaces. 
Since the fixed-node energy is an upper bound to the 
true ground state, these parameters can be optimized to 
give the best approximation to the ground-state wave 
function. In order to optimize these variational param- 
eters, we include them as slowly diffusing coordinates 
in a preliminary GFMC calculation. The parameters 
evolve slowly in imaginary time, equilibrating around the 
lowest-energy state consistent with the chosen form of the 
trial wave function. [4l| 

The ground-state energy Eq can be obtained from: 

(*v|ff|*o) _ (^oli^o) (23) 



En = 



(*y|*o) 



(*o|*o) 



Expectation values of quantities that do not commute 
with the Hamiltonian can be calculated using a combi- 
nation of the mixed and variational estimate: 



(*o|5|* > ~ 2(* |S|4V) - (^v\S\^ v ) , 



(24) 
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FIG. 2: (color online) Normalized neutron-matter variation- 
ally optimized pairing function 4>(r) for fcya = —10 (solid 
lines) and /cfa = — 1 (dotted lines), for different directions 
in the periodic simulation volume (in terms of rising expanse 
they correspond to the 001, Oil, and 111 directions in the 
box). 



for the opposite-spin f(r) and by the corresponding equa- 
tion for the same-spin fp(r). Since the f(r) and fp{r) 
we get are nodclcss, they do not affect the final result 
apart from reducing the statistical error. The fixed-node 
approximation guarantees that the result we obtain for 
one set of pairing function parameters in Eq. (|27[) will 
be an upper bound to the true grou nd-state energy of 
the system. As in previous works |32l . lilj the parameters 
are optimized in the full QMC calculation, providing the 
best possible nodal surface, in the sense of lowest fixed- 
node energy, with the given form of trial function. We 
utilize this upper-bound property to get as close to the 
true ground-state energy as possible. 

Given the finite-size analysis shown in Refs. [HI, l38j . 
we have performed all calculations for 66-68 particles in 
periodic boundary conditions. The equation of state is 
determined from the 66 particle results, and the pairing 
gap from the odd-even staggering. We have separately 
optimized the wave-function parameters at each density, 
and show the results for <p(r) (normalized each time to 
the value at zero separation) for the largest and smallest 
density we have considered (hp a = —10 and kpa = —1) 
in Fig. H 



where S is the operator corresponding to the relevant 
physical quantity, and the error in this expression is of 
second order in ^q— \EV- Such a combination of estimates 
is often called the "extrapolated estimate" . 

C. Trial wave function 

We take the trial wave function to be of the Jastrow- 
BCS form with fixed particle number: 

*y =I]>(*j) II fr-->v/i\[f<r,y:A \[ 0(r«O] 

i^J i'^f hi' '<■<]' 

(25) 

and periodic boundary conditions. The primed (un- 
primed) indices correspond to spin-up (spin-down) neu- 
trons. The pairing function 4>{r) is a sum over the mo- 
menta compatible with the periodic boundary conditions. 
In the BCS theory the pairing function is: 

and here it is parametrized with a short- and long-range 
part as in Ref. (ill ]: 

4>(r)=p(r)+ (27) 

n, I<I C 

where I = n 2 + n 2 y + n 2 z using the parameters defined in 
Eq. (fT2|) . The Jastrow part is taken from a lowest-order- 
constrained- variational method [42| calculation described 
by a Schrodinger-likc equation: 

--V 2 f(r)+v(r)f(r)=Xf(r) (28) 



D. Equation of state 

We first examine the equation of state of low-density 
neutron matter, in particular its evolution from the weak- 
to strong-coupling regime and the impact of adding p- 
wave interactions between the neutrons. In Fig. [3] we 
show the equation of state for low-density neutron matter 
with the s-wave interaction potential along with the new 
AV4 results. It is clear that when the density is very 
low, the s-wave contribution is dominant, and our results 
for the lowest densities remain unchanged. At higher 
densities the energy is higher with the contribution of 
the p-wave interaction. For the highest density examined, 
kpa = —10, this change is approximately 7%, while for 
kpa = —5 it is only 1%. Nonpcrturbative corrections at 
the highest density considered could reduce the difference 
between the s-wave interaction and AV4 results slightly. 
The curve at lower densities shows the analytic result 
[34[ described in section Hi Al Our calculations extend to 
lower densities than other microscopic calculations, and 
agree with the trend implied by the Lee- Yang result. 

We also compare our QMC AV4 results from Fig. [3] for 
the ground-state energy to other calculations extending 
to larger Fermi momenta. In Fig. 0] we compare our re- 
sults to (approximate) variational hypernetted-chain cal- 
culations by Friedman and Pandharipande Q, and an- 
other calculation by Akmal, Pandharipande, and Raven- 
hall (APR) @. We also include two Green's Function 
Monte Carlo results for 14 neutrons with more complete 
Hamiltonians (Toj . a result following from a Brueckner- 
Bethe-Goldstone expansion flSjf. a difermion EFT result 
(shown arc the error bands) [TlJ, the latest Auxiliary- 
Field Diffusion Monte Carlo (AFDMC) result (discussed 
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FIG. 3: (color online) Equation of state for neutron matter 
using different potentials. Shown are QMC results for the 
s-wave potential (circles) and for the AV4 (squares). Also 
shown is the analytic expansion of the ground-state energy of 
a normal fluid (line). 
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FIG. 4: (color online) Equation of state for neutron matter 
compared to various previous results. Despite quantitative 
discrepancies, all calculations give essentially similar results. 
Our lowest density corresponds to fen = — 1. 



below) (28|, a Dirac-Brueckner-Hartree-Fock calculation 
[l2j . a lattice chiral EFT method at next to leading or- 
der [lij (see also Ref . [l5[ ) , and an approach that makes 
use of chiral N 2 LO three- nucleon forces. [l6[ Of these, 
Refs. [§], [iU, and [l6[ include a three- nucleon inter- 
action, though at the densities we consider, these are not 
expected to be significant. Qualitatively all of these re- 
sults agree within 20%. 

A series of ab initio calculations for neutron matter us- 
ing the AFDMC method have been published beginning 
in 2005. [2^| After our analysis of the finite-size effects - 
described for BCS in section III Bl and for QMC in Refs. 
[H Hi| - was published in late 2007, the AFDM(^roup 
repeated their calculations for larger systems, [28|, |3(| 
bringing them closer to our results, though still, as can 
be seen from Fig. @] the results are distinct. Given the 
ab initio nature of the powerful AFDMC method, [43| we 
have attempted to compare results more extensively. The 
advantage of the AFDMC approach is that it includes an 
interaction which is more complete than the simpler ones 
used here. The disadvantage of the AFDMC approach is 
that it does not provide a variational bound to the energy, 
and hence the wave functions are chosen from another 
approach. In the calculations of Refs. [25|, HH, [3(| the 
wave function was taken from a Correlated-Basis Func- 
tion (CBF) approach that included a BCS-like initial 
state. The pairing in that variational state is unusually 
large, and in fact increases as a fraction of Ep when the 
density is lowered. 

The QMC AV4 results use a wave function that has 
been variationally optimized. QMC thus gives ener- 
gies that are considerably lower than the AFDMC re- 



sults. As both the wave functions and the interactions 
are different in the previous QMC and AFDMC results, 
we have repeated our calculations using the same input 
wave function [44| used by the AFDMC group (which 
comes from the same Correlated-Basis Function calcula- 
tion) at kF = 0.4 fm" 1 and at fcpa = —10. We find that 
in QMC the AV4 results for the optimized wave func- 
tion [0.5866(6) MeV and 0.5870(3) MeV, respectively] are 
consistently lower in energy than those using the CBF as 
input [0.6254(9) MeV and 0.6014(7) MeV, respectively]. 
This means that they are closer to the true ground-state 
energy for the Hamiltonian we consider. It would be 
worth studying in more detail the differences arising from 
the different Hamiltonians; the most important remain- 
ing differences are likely the spin-orbit and pion-exchange 
terms in the p-wave interaction. Extensions of previous 
GFMC calculations [lj| to lower densities would help to 
resolve these issues. 

It is interesting to note that at the lowest densities con- 
sidered, the AFDMC and QMC results are still distinct. 
At those densities contributions of p- and higher partial 
waves in the Hamiltonian should be very small, and thus 
the two methods should give identical results. The three- 
nucleon interaction included in the AFDMC calculations 
is one possible source of the difference, though this ap- 
pears unlikely at the smallest densities considered. This 
suggests that the CBF wave function at very low densi- 
ties is problematic; additional studies with Jastrow-BCS 
or other wave functions would be useful. 
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E. Pairing gap and quasiparticle spectrum 

We have also performed calculations for the zero- 
temperature pairing gap using the AV4 interaction. 
These follow from our knowledge of the ground-state en- 
ergy, through the use of the the odd-even staggering for- 
mula: 

A = E(N + 1) - i [E(N) + E(N + 2)] , (29) 

where N is an even number of particles. The results for 
the gap are shown in Fig. [5] The main conclusion is that 
the gap remains essentially unchanged with the inclusion 
of the p-wave interactions. Even at the highest density 
examined, kpa = —10, the gap is within statistical errors 
the same comparing s-wave and AV4 interactions. This 
implies that the dominant contributions to the gap come 
from the s-wave part of the interaction. 

Our results indicate that the gap is suppressed by 
approximately a factor of two from the BCS value at 
kpa = —1, roughly consistent with the Gorkov and 
Mclik-Barkhudarov, Eq. ([3]), polarization suppression. 
In cold atoms, this suppression from BCS is reduced as 
the density increases, with a smoothly growing fraction 
of the BCS results as wc move from the BCS to the 
BEC regime. At unitarity the measured pairing gaps 
[45|-l47j are 0.45(0.05) of the Fermi energy, for a ratio 
A/Abcs ~ 0.65,m agreement with predictions by QMC 
methods. [32T. KlL lijjj In neutron matter, though, the finite 
range of the potential reduces A/Ep as the density in- 
creases. We find a ratio A/Abcs that increases slightly 
from \kpa\ = 1 to 2.5, but thereafter remains roughly 
constant. 



TABLE I: Gap differences at various kpa calculated in per- 
turbation theory. Perturbative estimates based on AV4 cal- 
culations. 



kpa 


k F [fm" 1 ] 


A(AV4) [MeV] A(AV4) - A(s) [MeV] 


-5.0 


0.27 


0.48 (0.04) 0.012 (0.008) 


-7.5 


0.40 


0.77 (0.08) 0.11 (0.03) 


-10.0 


0.54 


1.05 (0.11) 0.16 (0.06) 



We also used our AV4 calculations to compute the dif- 
ference between s-wave and AV4 interaction gaps in per- 
turbation theory, in an attempt to isolate the effects of 
the addition of the p-wave interaction. This perturba- 
tion theory may not be accurate for the highest density 
considered, since the s-wave and AV4 ground states are 
somewhat different in energy. It should give an accurate 
picture at lower densities, though, and in particular iso- 
late the sign of the change arising from the p-wave terms 
in the interaction. Using perturbation theory yields much 
smaller statistical errors than comparing the separate s- 
wave and AV4 calculations. Table I shows that the p- 
wave interactions increase the pairing gap modestly over 
the range of densities considered. The p-wave interac- 
tions apparently decrease the magnitude of the polariza- 
tion corrections, though the change is only approximately 
15 % at the highest density considered. 

In Fig. [H] we compare our results to selected previous 
results: a Correlated-Basis Function calculation by Chen 
et al. [20| . an extension of the polarization-potential 
model by Wambach et al. 21[, a medium-polarization 
calculation by Schulzc et al. [22 j_a renormalization group 
calculation by Schwenk et al. 23J, a Brueckncr calcula- 
tion by Cao et al. [26| , a determinantal lattice QMC ap- 
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FIG. 5: (color online) Superfluid pairing gap versus feo for 
neutron matter using different potentials. Shown are QMC 
results for the s-wave potential (circles) and for the AV4 
(squares). Also shown is the mean-field BCS result (line). 
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FIG. 6: (color online) Superfluid pairing gap versus kpa for 
neutron matter compared to previous results. 
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FIG. 7: (color online) The neutron-matter energy of the quasi- 
particle excitations of the system in QMC AV4 (squares) ver- 
sus (k/kp) 2 at kpa = —10. Also shown are the BCS contin- 
uum results (line) as well as the QMC quasiparticle spectrum 
that follows from a simple s-wave Hamiltonian (circles). 



proach [29( , and finally the newer Correlated-Basis Func- 
tion calculation by Fabrocini et al. [2o| that was used as 
an input wave function in the two AFDMC calculations 
of 2005 and 2008. [H, H 

The results of our calculations are much larger than 
most diagrammatic [20l - l22j and renormalization group 
[23j approaches. As these approaches assume a well- 
defined Fermi surface or calculate polarization correc- 
tions based on single-particle excitations it is not clear 
how well they can describe neutron matter in the strongly 
paired regime, or the similar pairing found in cold atoms. 
Ref. prjj . which incorporates self-energy corrections and 
screening at the RPA level within Brucckner theory, ap- 
pears to give results similar to ours. However, these val- 
ues disagree with our lower-density results and, perhaps 
more importantly, at the lowest density reported the gap 
is larger than the mean- field BCS value (see section Ul Al) . 
On a similar note, Refs. [2l[ and [23[ make use of a weak- 
coupling formula to calculate the pairing gap, similarly 
to the Eqs. and we discussed in section iHAl The 
prefactor they use is justified based on predictions in the 
theory of 3 He. However, the concept itself of a Fermi 
surface is not well-defined in these strongly paired sys- 
tems: in 3 He, in contrast to the present case, the gap is 
considerably smaller than the Fermi energy. 

Our results are also somewhat different from the 
AFDMC results of Ref. (28|. We have once again re- 
peated our QMC calculations for the gap using the CBF 
wave function as input. We find something quite interest- 
ing: the QMC method using the AV4 potential and the 
CBF input wave function at kp = 0.4 fm _1 (which gave 
an energy higher than the variationally optimized input 
wave function, see section H*IID[) gives a gap of 1.21(17) 
McV, thus reproducing the AFDMC result, which uses 



the same input wave function and the much more com- 
plicated AV8'+UIX interaction, this being 1.45(15) MeV. 
This too suggests that the most important contributions 
to the gap come from the s-wave part of the interac- 
tion. On the other hand, our results seem to qualitatively 
agree (at least for the lowest densities considered) with a 
determinantal Quantum Monte Carlo lattice calculation 
[29l | which, however, only includes the s-wave component 
in the interaction. This may imply that a consensus is 
emerging, in that both these calculations find a gap that 
is suppressed with respect to the mean-field BCS result 
but is still a substantial fraction of the Fermi energy. 
Finally, let us note that, as mentioned before when dis- 
cussing Fig. [SJ the AV4 results for the optimized wave 
function are very similar to those using an s-wave poten- 
tial. 

We have also calculated the quasiparticle excitation 
spectrum using the AV4 interaction. The minimum of 
these results provides the pairing gap in Fig. [5] In Fig. 
[7]we show both the s-wave Hamiltonian results, as well as 
the AV4 results. In cold atoms at unitarity and beyond 
(the BEC regime) the quasiparticle minimum energy is 
at a momentum significantly smaller than the Fermi 
momentum. Here, though, the minimum corresponds 
closely to the neutron Fermi momentum. Although the 
QMC minimum (pairing gap) is much smaller than the 
BCS minimum, the dispersion around the minima is quite 
similar. Just like in the case of cold atoms, microscopic 
results for the quasiparticle energy spectra [48| can be 
used to constrain density functional calculations. [49| 



F. Distribution functions 

Using the QMC AV4 interaction, we have also cal- 
culated distribution functions. In Fig. [5] we show the 
momentum distribution at three densities, calculated as 
the Fourier transform of the one-body density matrix, 
through: 



*y(ri,...,r„) 



, (30) 



where the curly brackets denote a stochastic integration 
over the angles and we perform the integral over Sr = 
r' n — r„ I on a line using Gaussian quadratures to avoid 
statistical errors due to the oscillatory radial dependence. 

As expected from standard BCS theory, we see that 
the spread of the momentum distribution around fi is 
approximately 2A. For kpa = —1 the momentum dis- 
tribution looks very similar to that of a free Fermi gas. 
At large \kpa\ (when the gap is approximately half the 
Fermi energy) this fact implies that there is no clearly 
defined Fermi surface. We note that in Fig. [5] the results 
for kpa = —5 seems to be more "broadened" than that 
of kpa = —10, even though as we can see in Fig. [6] the 
gap at kpa = —5 is considerably smaller than the one 
at kpa = —10. This is easily resolved if one looks at 
the pairing gap not in absolute units (MeV) but divided 
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FIG. 8: (color online) The neutron-matter momentum dis- 
tribution in QMC versus (k/kp) 2 at kpa = —1 (squares), 
kpa — — 5 (diamonds), and kpa = —10 (circles). Also shown 
are the continuum BCS results at kpa = —1 (dashed line), 
kpa = —5 (dotted line), and kpa = —10 (solid line). 
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FIG. 9: (color online) The neutron-matter pair-distribution 
function in QMC as a function of the distance times the Fermi 
momentum at kpa = —10. The distribution functions given 
are the g c {r) for opposite-spin (circles) and same-spin pairs 
(squares), as well as the gu {r) for opposite-spin pairs (trian- 
gles). 



with the Fermi energy, as shown in Fig. [5] The case of 
kpa = —5 has a bigger gap in units of the Fermi energy, 
and that is what leads to the observed behavior in the 
momentum distributions shown in Fig. [5] 

We have also computed the pair-distribution functions 
at kpa = — 10 using the AV4 potential, and have plotted 
them in Fig. [S] These are calculated from expectation 



g P (r)=AYm \6(r ij -r)O[ i \* v ) 



i<j 



(31) 



where we are initially interested in the case in which the 
operator is simply unity, and the normalization factor 
A is such that gi(r) = g c (r) goes to one at large dis- 
tances. These pair-distribution functions provide sum 
rules related to density- and other response functions 
versus density and momentum. The solid line in the fig- 
ure shows the pair-distribution function of noninteracting 
(NI) fermions with parallel spins: 



,ni 



(r) = 1 - 



{k F rf 



[s'm(kpr) — kpr cos(/ci?r)] 2 . (32) 



The noninteracting result is very close to the QMC sim- 
ulation for same-spin particles. On the other hand, the 
value of the opposite-spin distribution function is very 
small at short distances, reflecting the repulsive core in 
the AV18 potential. Since our interaction is more com- 
plicated than a simple s-wave component, Eq. (|31|) can 
also be applied to the Majorana exchange operator P M 
which was used in Eq. (|18[) . The distribution function 
for that operator is also shown in Fig. |9]for the density of 
interest. It tracks the behavior of the central (unit oper- 
ator) distribution function for short distances, but then 
reduces to approximately half the standard distribution 
at kpr « 1.7. 



IV. CONCLUSIONS 

To conclude, we have calculated the equation of state 
and pairing gap of low-density neutron matter with the 
AV4 interaction from kp = 0.054 to 0.54 fm _1 , corre- 
sponding to \kpa\ from 1 to 10. The calculated equation 
of state and pairing gap match smoothly with the known 
analytic results at low densities, and provide important 
constraints in the strong-coupling regime at large kpa. 
We have also calculated the quasiparticlc spectrum, mo- 
mentum distribution, and pair-distribution functions for 
low-density neutron matter. The low-density equation of 
state can help constrain Skyrme mean-field models of fi- 
nite nuclei. The pairing gap for low-density neutron mat- 
ter is relevant to Skyrme-Hartree-Fock-Bogoliubov cal- 
culations 0] of neutron-rich nuclei and to neutron-star 
physics, since it is expected to influence the behavior of 
the crust. @ 

More specifically, a magnetic field in the neutron star 
crust would have to be approximately 10 G to over- 
come this gap and thus polarize neutron matter; such a 
value of the magnetic field is not implausible within the 
context of magnetars. Similarly, the fact that the mag- 
nitude of the gap is not as small as previously expected 
implies that a new mechanism that makes use of super- 
fluid phonons is competitive to the heat conduction by 
electrons in magnetized neutron stars. @ 
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Our results for the gap at the low-density regime, fol- 
lowing from a variationally optimized approach that in- 
cludes the dominant well-known terms in the Hamilto- 
nian, can function as a benchmark with which other cal- 
culations can be compared. Equally important, future 
experimental tests in cold atoms, at least in the very low 
density regime up to |fcpa| = 2 appear to be within the 
reach of possibility. Similar comparisons may be made 
with other observables including the pair-distribution 
function and momentum distribution. We believe that 
these calculations of the equation of state, pairing gap, 
and quasiparticle dispersion can be used as constraints 
of nuclear density functionals. 
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